<?php

    /**
     * @package JAMA
     *    For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
     *    unit lower triangular matrix L, an n-by-n upper triangular matrix U,
     *    and a permutation vector piv of length m so that A(piv,:) = L*U.
     *    If m < n, then L is m-by-m and U is m-by-n.
     *    The LU decompostion with pivoting always exists, even if the matrix is
     *    singular, so the constructor will never fail. The primary use of the
     *    LU decomposition is in the solution of square systems of simultaneous
     *    linear equations. This will fail if isNonsingular() returns false.
     * @author  Paul Meagher
     * @author  Bartosz Matosiuk
     * @author  Michael Bommarito
     * @version 1.1
     * @license PHP v3.0
     */
    class PHPExcel_Shared_JAMA_LUDecomposition {
        const MATRIX_SINGULAR_EXCEPTION = "Can only perform operation on singular matrix.";
        const MATRIX_SQUARE_EXCEPTION   = "Mismatched Row dimension";
        /**
         *    Decomposition storage
         * @var array
         */
        private $LU = [];
        /**
         *    Row dimension.
         * @var int
         */
        private $m;
        /**
         *    Column dimension.
         * @var int
         */
        private $n;
        /**
         *    Pivot sign.
         * @var int
         */
        private $pivsign;
        /**
         *    Internal storage of pivot vector.
         * @var array
         */
        private $piv = [];

        /**
         *    LU Decomposition constructor.
         * @param $A Rectangular matrix
         * @return Structure to access L, U and piv.
         */
        public function __construct($A) {
            if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {
                // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
                $this->LU = $A->getArray();
                $this->m  = $A->getRowDimension();
                $this->n  = $A->getColumnDimension();
                for ($i = 0; $i < $this->m; ++$i) {
                    $this->piv[$i] = $i;
                }
                $this->pivsign = 1;
                $LUrowi        = $LUcolj = [];
                // Outer loop.
                for ($j = 0; $j < $this->n; ++$j) {
                    // Make a copy of the j-th column to localize references.
                    for ($i = 0; $i < $this->m; ++$i) {
                        $LUcolj[$i] = &$this->LU[$i][$j];
                    }
                    // Apply previous transformations.
                    for ($i = 0; $i < $this->m; ++$i) {
                        $LUrowi = $this->LU[$i];
                        // Most of the time is spent in the following dot product.
                        $kmax = min($i, $j);
                        $s    = 0.0;
                        for ($k = 0; $k < $kmax; ++$k) {
                            $s += $LUrowi[$k] * $LUcolj[$k];
                        }
                        $LUrowi[$j] = $LUcolj[$i] -= $s;
                    }
                    // Find pivot and exchange if necessary.
                    $p = $j;
                    for ($i = $j + 1; $i < $this->m; ++$i) {
                        if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
                            $p = $i;
                        }
                    }
                    if ($p != $j) {
                        for ($k = 0; $k < $this->n; ++$k) {
                            $t                = $this->LU[$p][$k];
                            $this->LU[$p][$k] = $this->LU[$j][$k];
                            $this->LU[$j][$k] = $t;
                        }
                        $k             = $this->piv[$p];
                        $this->piv[$p] = $this->piv[$j];
                        $this->piv[$j] = $k;
                        $this->pivsign = $this->pivsign * -1;
                    }
                    // Compute multipliers.
                    if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
                        for ($i = $j + 1; $i < $this->m; ++$i) {
                            $this->LU[$i][$j] /= $this->LU[$j][$j];
                        }
                    }
                }
            } else {
                throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ARGUMENT_TYPE_EXCEPTION);
            }
        }    //    function __construct()

        /**
         *    Get lower triangular factor.
         * @return array Lower triangular factor
         */
        public function getL() {
            for ($i = 0; $i < $this->m; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    if ($i > $j) {
                        $L[$i][$j] = $this->LU[$i][$j];
                    } elseif ($i == $j) {
                        $L[$i][$j] = 1.0;
                    } else {
                        $L[$i][$j] = 0.0;
                    }
                }
            }
            return new PHPExcel_Shared_JAMA_Matrix($L);
        }    //    function getL()

        /**
         *    Get upper triangular factor.
         * @return array Upper triangular factor
         */
        public function getU() {
            for ($i = 0; $i < $this->n; ++$i) {
                for ($j = 0; $j < $this->n; ++$j) {
                    if ($i <= $j) {
                        $U[$i][$j] = $this->LU[$i][$j];
                    } else {
                        $U[$i][$j] = 0.0;
                    }
                }
            }
            return new PHPExcel_Shared_JAMA_Matrix($U);
        }    //    function getU()

        /**
         *    Alias for getPivot
         * @see getPivot
         */
        public function getDoublePivot() {
            return $this->getPivot();
        }    //    function getPivot()

        /**
         *    Return pivot permutation vector.
         * @return array Pivot vector
         */
        public function getPivot() {
            return $this->piv;
        }    //    function getDoublePivot()

        /**
         *    Count determinants
         * @return array d matrix deterninat
         */
        public function det() {
            if ($this->m == $this->n) {
                $d = $this->pivsign;
                for ($j = 0; $j < $this->n; ++$j) {
                    $d *= $this->LU[$j][$j];
                }
                return $d;
            } else {
                throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MATRIX_DIMENSION_EXCEPTION);
            }
        }    //    function isNonsingular()

        /**
         *    Solve A*X = B
         * @param  $B  A Matrix with as many rows as A and any number of columns.
         * @return  X so that L*U*X = B(piv,:)
         * @PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.
         * @PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.
         */
        public function solve($B) {
            if ($B->getRowDimension() == $this->m) {
                if ($this->isNonsingular()) {
                    // Copy right hand side with pivoting
                    $nx = $B->getColumnDimension();
                    $X  = $B->getMatrix($this->piv, 0, $nx - 1);
                    // Solve L*Y = B(piv,:)
                    for ($k = 0; $k < $this->n; ++$k) {
                        for ($i = $k + 1; $i < $this->n; ++$i) {
                            for ($j = 0; $j < $nx; ++$j) {
                                $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                            }
                        }
                    }
                    // Solve U*X = Y;
                    for ($k = $this->n - 1; $k >= 0; --$k) {
                        for ($j = 0; $j < $nx; ++$j) {
                            $X->A[$k][$j] /= $this->LU[$k][$k];
                        }
                        for ($i = 0; $i < $k; ++$i) {
                            for ($j = 0; $j < $nx; ++$j) {
                                $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
                            }
                        }
                    }
                    return $X;
                } else {
                    throw new PHPExcel_Calculation_Exception(self::MATRIX_SINGULAR_EXCEPTION);
                }
            } else {
                throw new PHPExcel_Calculation_Exception(self::MATRIX_SQUARE_EXCEPTION);
            }
        }    //    function det()

        /**
         *    Is the matrix nonsingular?
         * @return true if U, and hence A, is nonsingular.
         */
        public function isNonsingular() {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($this->LU[$j][$j] == 0) {
                    return false;
                }
            }
            return true;
        }
    }
